Install libraries if there are no already in your computer.
## [1] "raster"
## Loading required package: raster
## Warning: package 'raster' was built under R version 4.0.5
## Loading required package: sp
## [1] "terra"
## Loading required package: terra
## Warning: package 'terra' was built under R version 4.0.5
## terra version 1.1.17
## [1] "snow"
## Loading required package: snow
## Warning: package 'snow' was built under R version 4.0.4
## [1] "parallel"
## Loading required package: parallel
##
## Attaching package: 'parallel'
## The following objects are masked from 'package:snow':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, clusterSplit, makeCluster, parApply,
## parCapply, parLapply, parRapply, parSapply, splitIndices,
## stopCluster
## [1] "foreach"
## Loading required package: foreach
## Warning: package 'foreach' was built under R version 4.0.5
## [1] "rasterVis"
## Loading required package: rasterVis
## Warning: package 'rasterVis' was built under R version 4.0.5
## Loading required package: lattice
## Loading required package: latticeExtra
## Warning: package 'latticeExtra' was built under R version 4.0.5
## [1] "mapview"
## Loading required package: mapview
## Warning: package 'mapview' was built under R version 4.0.4
## GDAL version >= 3.1.0 | setting mapviewOptions(fgb = TRUE)
## [1] "ggplot2"
## Loading required package: ggplot2
##
## Attaching package: 'ggplot2'
## The following object is masked from 'package:latticeExtra':
##
## layer
## [1] "lubridate"
## Loading required package: lubridate
##
## Attaching package: 'lubridate'
## The following objects are masked from 'package:terra':
##
## intersect, union
## The following objects are masked from 'package:raster':
##
## intersect, union
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
## [1] "pryr"
## Loading required package: pryr
## Warning: package 'pryr' was built under R version 4.0.5
## Registered S3 method overwritten by 'pryr':
## method from
## print.bytes Rcpp
##
## Attaching package: 'pryr'
## The following object is masked from 'package:terra':
##
## dots
## The following object is masked from 'package:raster':
##
## subs
Download files, alternative is possible use the function raster::getData. More details on https://worldclim.org/data/worldclim21.html
#download worldclim climate data
download.file('https://biogeo.ucdavis.edu/data/worldclim/v2.1/base/wc2.1_10m_tavg.zip',temp, mode="wb")
#raster::getData()
unzip(temp, exdir=tempd)
#dir(tempd) # explore the temp directory
Explore the data
file_names = list.files(tempd,'wc2.1_10m_tavg')
tmean <- stack(file.path(tempd, file_names))
mapview(tmean[[7]])
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ;
## the supplied Raster* has 2332800
## ... decreasing Raster* resolution to 5e+05 pixels
## to view full resolution set 'maxpixels = 2332800 '
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition
## Warning in rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, : 104
## projected point(s) not finite
## Warning in rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, : 104
## projected point(s) not finite
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition
## Paraleling Create cluster and load functions. Alternative you can create your own scrip file with function and just load this file.
cl <- makeCluster(4, type = "SOCK")
out=clusterEvalQ(cl,library(raster))
clusterExport(cl, list('my_function'))
#out =clusterEvalQ(cl, source('my_functions.R'))
We going to use two way of paralleling in R with the snow/parallel and foreach libraries.
system.time({
clusterExport(cl,list('tmean'))
mean_planet = parLapply(cl,c(1:12), function(l){
my_function(tmean[[l]])
})
mean_planet = unlist(mean_planet)
})
## user system elapsed
## 0.02 0.01 2.05
system.time({
mean_planet2 = foreach(l = c(1:12),
.combine = rbind,
.errorhandling = 'remove') %dopar%
my_function(tmean[[l]])
})
## Warning: executing %dopar% sequentially: no parallel backend registered
## user system elapsed
## 3.48 0.53 4.02
results = data.frame(layer_name = names(tmean),
month = month(1:12,label=TRUE, abbr = TRUE),
parLapply_out =mean_planet,
foreach_out = mean_planet2)
The first plot show the global mean temperature for each month. The second figure show that the two paralleling outcome the same results
ggplot(results,aes(x=month,y=parLapply_out, col='parLapply_out'))+
geom_point()+labs(y='Mean temperature (°C)', col='', x='')+
geom_point(aes(y=foreach_out,col = 'foreach_out'))+theme_minimal()
We download the raster and make a split overlapped extent. The overlapped extent is only for the use of focus function, any other process should use a non-overlapped extent! Also is possible use the function focus in raster. but it is slower function.
file_raster = 'https://www.dropbox.com/s/onp3f75pab7z3j9/intro-1508529851.jpg?dl=1'
r_raster = stack(file_raster)
e =extent(r_raster)
extent_r =list(r1 =extent(0, floor(e@xmax/2)+2, 0, floor(e@ymax/2)+2),
r2 = extent(ceiling (e@xmax/2), e@xmax, floor (e@ymax/2), e@ymax),
r3 = extent(0, floor (e@xmax/2)+2, floor (e@ymax/2), e@ymax),
r4 = extent(ceiling (e@xmax/2), e@xmax, 0, floor (e@ymax/2)+2))
We call the parLapply from a lapply environment, so we need export elements from this environment (env).
tmean_smooth = lapply(c(1:3) , function(l){
env= where('l')
clusterExport(cl, list('extent_r','r_raster','l'),envir=env)
r_parts = parLapply(cl,seq_along(extent_r),function(part){
r = r_raster[[l]]
r_part <- crop(r, extent_r[[part]])
#focal(r_part,w=matrix(1,3,3), fun='mean',na.rm=TRUE)
terra::focal(r_part,w=matrix(1,3,3), fun='mean',na.rm=TRUE)
#r_part
})
m1 <- mosaic(r_parts[[1]], r_parts[[2]],r_parts[[3]],r_parts[[4]], fun=mean)
m1
})
#r_smooth = stack(tmean_smooth[[1]],tmean_smooth[[2]],tmean_smooth[[3]])
r_smooth = do.call(stack,tmean_smooth)
plotRGB(r_smooth)
stopCluster(cl)
Versus original
plotRGB(r_raster)